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We lay the foundations for the construction of analytic expressions for Fourier-domain gravita- 
tional waveforms produced by eccentric, inspiraling compact binaries in a post-circular or small- 
eccentricity approximation. The time-dependent, "plus" and "cross" polarizations are expanded 
in Bessel functions, which are then self-consistently re-expanded in a power series about zero ini- 
tial eccentricity to eighth order. The stationary phase approximation is then employed to obtain 
explicit analytic expressions for the Fourier transform of the post-circular expanded, time-domain 
signal. We exemplify this framework by considering Newtonian-accurate waveforms, which in the 
post-circular scheme give rise to higher harmonics of the orbital phase and amplitude corrections 
both to the amplitude and the phase of the Fourier domain waveform. Such higher harmonics lead 
to an effective increase in the inspiral mass reach of a detector as a function of the binary's eccen- 
tricity Co at the time when the binary enters the detector sensitivity band. Using the largest initial 
eccentricity allowed by our approximations (eo < 0.4), the mass reach is found to be enhanced up 
to factors of approximately 5 relative to that of circular binaries for Advanced LIGO, LISA, and 
the proposed Einstein Telescope at a signal-to-noise ratio of ten. A post-Newtonian generalization 
of the post-circular scheme is also discussed, which holds the promise to provide "ready-to-use" 
Fourier-domain waveforms for data analysis of eccentric inspirals. 

PACS numbers: 04.30.-w,04.30.Db,04.30.Tv,04.25.Nx 



I. INTRODUCTION 

The detection and characterization of gravitational 
waves (GWs) hold the promise to reveal previously 
unattainable, yet very valuable astrophysical informa- 
tion (sec [l[ for a recent review). Ground-based de- 
tectors, such as the Laser Interferometer Gravitational 
Wave Observatory (LIGO) 0, VIRGO [|, GEO [|, and 
TAMA d, have started acquiring data at or near design 
sensitivity. The space-borne Laser Interferometer Space 
Antenna (LISA) [6| may be launched within the next 
decade, while future Earth-based third-generation detec- 
tors, such as the proposed Einstein Telescope (ET) Q, 
are currently being planned. These detectors are ex- 
pected to observe several different astrophysical GW 
sources, one of the most promising of which arc compact 
binary inspirals. 

Black hole (BH) binaries are considered one of the 
main targets for GW detection and their evolution can 
be roughly divided into an inspiral phase and a merger 
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plus ringdown phase. Peters and Mathews showed that 
eccentric inspirals circularize via GW emission [1, 0], 
and thus, it was traditionally thought that eccentricity 
would not play a major role in GW detection and anal- 
ysis. One can show, for example, that to leading order 
e/eo (///o)~^^/^^ [3, which implies that if a binary 
enters the sensitivity band of a ground-based interfer- 
ometer at 20 Hz with an initial eccentricity of 0.1, its 
eccentricity is reduced to 0.01 before the system reaches 
200 Hz. Circularization via GW emission, however, is not 
absolute, since systems that enter a detector's sensitiv- 
ity band with large enough eccentricity can retain some 
residual eccentricity before they merge or exit the band. 
For example, if a binary enters the sensitivity band at 
20 Hz with an initial eccentricity of 0.4, its eccentricity 
remains significant while in band, and is reduced to 0.01 
only by the time the frequency reaches 10^ Hz. 

Astrophysical scenarios have been proposed that pre- 
dict that binary inspiral signals could enter the sensitivity 
band of GW detectors with non-negligible eccentricity. 
More precisely. Earth-based detectors are expected to be 
sensitive to stellar mass BH/BH binaries, which might 
fail to completely circularize before merger, leading to po- 
tentially large eccentricities in the detector band [lll.[l3|. 
Wen [ll| studied the evolution of the inner binary of 
triple BH systems in globular clusters. She found that 
approximately 30% of these binaries could merge via the 
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Kozai mechanism with eccentricities > 0.1 by the time 
they enter the LIGO band at frequency ~ 10 Hz. O'Leary 
et al. studied the density cusp of stellar mass BHs 
that forms around supermassive BHs in galactic nuclei 
because of mass segregation. They found that, in such 
dense environments, hyperbolic BH-BH encounters can 
lead to the formation of bound systems, and that most 
of these binaries (~ 90%) would have eccentricity > 0.9 
when they enter the LIGO band. Several classes of com- 
pact binary sources for LISA are also predicted to be 
eccentric when they enter the detector's band. In or- 
der to improve readability of this paper, we have rele- 
gated a brief review of these astrophysical scenarios to 
Appendix [A] 

For binaries with non-negligible eccentricities, a dedi- 
cated search using matched filtering techniques and ec- 
centric orbit templates would be necessary for detection 
and extraction of astrophysical information. Until now, 
however, closed-form analytic expressions for such tem- 
plates have been lacking, with most studies concentrating 
on the circular case. The modeling of eccentric orbits is 
much more difficult, since it requires knowledge not only 
of the orbital phase and frequency, but also of frequen- 
cies associated with the higher harmonics of the eccentric 
motion, as well as of frequencies associated with PN pre- 
cession effects, such as pericenter precession. 

Peters and Mathews' seminal work [l3|, and Peters' 
follow-up calculation of the angular momentum flux and 
evolution of orbital elements at Newtonian order [T^ . 
laid the foundations for the calculation of energy and an- 
gular momentum fluxes at higher PN orders [3, [l^, [13, 
Ej 0j [13, HH, J the latter mostly focusing on circular 
orbits. Analytic GW templates for circular binaries are 
now available at 3.5PN order in the pha se il,[ll,[2^ and 
3PN order in the amplitude (2^, [27l. l28l. [29 1 , and their as- 
sociated Fourier transforms have been computed in the 
stationary phase approximation. Recently Ref. [s^l pro- 
vided a method to construct high accuracy templates for 
elliptical binaries in the time domain by explicitly com- 
puting the post-adiabatic short period contributions at 
2.5PN order, to be added to the post-Newtonian expres- 
sions for the GW polarizations. A 3.5PN generalization 
of these templates was discussed in Ref. [31[ . 

Relatively few investigations of data analysis issues for 
eccentric inspirals have been performed. Martel and Pois- 
son [s^l quantified the accuracy to which circular orbit 
templates could capture signals from eccentric binaries 
in the LIGO band, finding that the signal-to-noise ra- 
tio (SNR) loss is significant for eccentricities above 0.1. 
Seto [33j studied parameter estimation in the context of 
eccentric gala ctic neutron star (NS) binaries, and Be- 
nacquista |34l. [35| carried out the first statistical inves- 
tigation of the harmonic structure of eccentric binary 
waveforms and their relevance for LISA GW detection. 
The analysis by Martel and Poisson was recently re- 
visited [la, [131, emphasizing the need for eccentric bi- 
nary templates. All these investigations, however, con- 
centrated on either time-domain waveforms or numerical 



Fourier transforms of such waveforms, which might not 
be desirable for data analysis purposes. 

The aim of this paper is to lay the foundations of a 
post-circular approximation that allows for the construc- 
tion of analytic, "ready-to-use" Fourier-domain wave- 
forms for eccentric binary inspirals. This approximation 
consists of expanding time-domain gravitational wave- 
forms in the eccentricity parameter eo, which is defined 
to be the eccentricity when the GW signal enters the de- 
tector sensitivity band. The resulting expression in then 
Fourier transformed in the stationary-phase approxima- 
tion. This scheme is an extension of the program initiated 
by Krolak, Kokkotas and Schafer several years ago [11| 
and it is intended to supplement the PN approximation, 
thus yielding a double or bivariate expansion in both the 
velocity of the binary members and the initial eccentric- 
ity. Although the PN scheme does not require the post- 
circular approximation for the construction of numerical 
Fourier-domain templates, closed-form analytic expres- 
sions for these templates can only be obtained through 
the incorporation of the post-circular approximation. 

The usefulness of "ready-to-use," analytic, frequency- 
domain waveforms is two-fold. Analytic expressions al- 
low us to study the structure of the eccentricity induced 
corrections to the Fourier transform of the signal; in 
turn, this structure allows us to explain features in the 
SNR that would otherwise be hidden by numerics. Sec- 
ondly, analytic expressions allow for fast implementa- 
tions of dedicated matched-filtering searches in a data 
analysis algorithm and allow us to sidestep fast Fourier 
transforms, which would drain numerical resources from 
Fourier-domain data analysis pipelines. 

We exemplify the post-circular approximation by con- 
sidering Newtonian expressions for the two GW polar- 
ization states of elliptic binaries. The cosines and sines 
of the GW phase are expanded in a truncated Bessel se- 
ries, whose argument is proportional to the eccentricity 
parameter. We find that the first 9 terms in the sum 
suffice to approximate the phase to better than 0.1% for 
eccentricities smaller than 0.4 (see Sec. [H] below). The 
Bessel series is then re-expanded to eighth order in e <C 1, 
which we find sufficient to capture the essential features 
of orbital dynamics for such small eccentricities. 

The structure of the time-domain gravitational wave- 
form for eccentric inspirals takes the form 

10 

e=i 

where F is the orbital frequency, I is the mean anomaly 
and C-f.x and S'+,x ai'e eighth-order power series in the 
eccentricity [see Eqs. ([X7|) - ([XTU)) and Eqs. ([BT|) - (|B36)) ]. 
The eccentricity is itself a function of the orbital fre- 
quency, which we invert in the limit e <C 1 and insert 
into the time-domain waveforms to obtain explicit ex- 
pressions whose only independent variables are the or- 
bital frequency and the mean anomaly. The prescription 
described above to compute time-domain waveforms is 
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a completion of the analytic work of Moreno-Garrido, 
Buitrago and Mediavilla [H, , improved through the 
resummation methods of Pierro and Pinto [4l|, |42| and 
the more general waveform expressions of Martel and 
Poisson [32j . 

Once closed-form, analytic expressions for the time- 
domain waveforms are obtained, we compute their 
Fourier transform in the stationary phase approxima- 
tion (SPA). This approximation derives from the asymp- 
totic method of integration by steepest descent, which 
allows one to systematically include higher harmonics in 
the frequency-domain waveforms. This higher-harmonic 
structure is found to fit perfectly in the formalism of 
Refs. [i^lSl, which was developed to account for higher 
harmonics due to PN amplitude corrections in circular- 
orbit binary waveforms. 

The Fourier transform of the response function is then 
found to take the form 



h^Af 
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(1.2) 



where A is an overall amplitude that depends on the 
system parameters (such as the masses of the binary 
members), while / is the dominant (quadrupole) GW fre- 
quency. The amplitudes and the phases ipe are small- 
eccentricity expansions [see Eqs. (|4.28p and (|C1|) ]. The 
^f's depend on the antenna pattern functions J^+.x [IHIj 
the initial eccentricity Cq and the GW frequency. The 
expansion of the phase reads 
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where /o is the frequency at which the eccentricity equals 
eo, which we choose to coincide with the low- frequency 
cut-off of the detector sensitivity band. This is in agree- 
ment with the Newtonian limit of Eq. (A9) of [s^ up to 
Oiel). 

One of the benefits of obtaining closed-form, analytic 
expressions for the Fourier transform of the waveforms is 
that its harmonic structure and its eccentricity-induced 
amplitude corrections become explicit. We find that 
these higher-harmonic eccentricity corrections increase 
the SNR for large total masses, in analogy to what was 
found for PN amplitude-corrected circular orbit wave- 
forms in [H, m, Figure [T] compares the optimal 
SNR of equal-mass binaries with eccentricity ep = and 
eo = 0.3 at the initial frequency of the sensitivity band 
of Advanced LIGO (AdvLIGO), FT and LISA (20 Hz, 
1 Hz or 10 Hz, and 10"'* Hz respectively). The source 
is located at distances of 100 Mpc for AdvLIGO and 
FT, and 3 Gpc for LISA. By "optimal" we mean the 
SNR measured by an observer located in a direction per- 
pendicular to the orbital plane (more precisely, we set 
t = /3 = 6's = (/)s = "05 = in the notation of Sect ion HIT)) . 
This SNR increase is rather generic, irrespective of the 
location of the source in the sky. 
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FIG. 1: SNR for an equal-mass binary at optimal orientation 
as a function of total mass for circular binaries and elliptic 
binaries with initial eccentricity of 0.3. The assumed low- 
frequency cut-offs for AdvLIGO, ET and LISA are 20 Hz, 
either 1 or 10 Hz, and 10~* Hz, respectively. The sources are 
at 100 Mpc for AdvLIGO and ET and at 3 Gpc for LISA. The 
initial eccentricity corresponds to that at the low-frequency 
cut-off. 
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FIG. 2: Normalized mass reach enhancement, as a function 
of initial eccentricity. The normalization is given by the value 
of the mass reach for circular binaries, namely Mo = 219M0, 
Mo = 44OM0(4397Mo) and Mo = 4.239 x lO^M© for Ad- 
vLIGO, ET and LISA, respectively. 



The inclusion of eccentricity in the waveforms leads to 
an increase in the mass reach as compared to circular 
waveforms. This is shown in Fig. [21 where we plot the 
mass reach enhancement M(eo)/Mo, where Mq = M(0). 
The mass reach M(eo) is here defined as the mass corre- 
sponding to an optimal SNR of ten, roughly correspond- 
ing to the largest mass visible to the detector. The mass 
reach usually increases with Cq, up to factors of order five 
for binaries with cq — 0.4. This result should still hold 
when PN corrections are included, since we expect these 
corrections to increase the mass reach. We conclude then 
that LISA could potentially observe moderately eccentric 
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binaries with total masses on the order of lO^M©. These 
results are in agreement with preliminary results from 
numerical relativity simulations of merging eccentric bi- 
naries [131 (see also (isj). 

An increase in mass reach in turn implies that, for a 
system of fixed total mass, the distance to which the sys- 
tem can be observed also increases with eccentricity. For 
example, IMBH mergers of total mass ~ 200Mo with or- 
bital eccentricity of eo — 0.3 would be observable by Ad- 
vLIGO up to approximately 1.26 Gpc (z ~ 0.26) with an 
SNR of 10, while SMBHs of total mass ~ 4 x lO^Mg and 
initial eccentricity cq — 0.3 would be visible by LISA up 
to approximately 100 Gpc (z ~ 10) with an SNR of 100. 
Such an increase in distance corresponds to an increase 
in accessible volume of up to a factor of approximately 
one hundred for systems with eo — 0.4. 

The post-circular approximation also allows us to de- 
termine how many harmonics are needed to reproduce 
the Fourier transform of the eccentric signal to some ac- 
curacy. For a system with eo = 0.01, we find that keeping 
up to the second or third harmonic suffices to reproduce 
the SNR of a signal that includes ten harmonics to 0(1) 
and O(10~^) respectively in the entire mass range. For a 
system with eo = 0.1, however, we find that a compara- 
ble accuracy in SNR requires including up to the fourth 
and fifth harmonic, respectively. Such an analysis allows 
us to conclude that the inclusion of up to the fourth har- 
monic suffices for SNR calculations when cq < 0.1, while 
for systems with 0.1 < eo < 0.3, one must really include 
eight harmonics or more. 

Although the waveforms we consider in our post- 
circular scheme are not accurate enough for a rigorous 
data analysis study, as we have ignored FN effects, they 
do provide insight as to the effect of eccentricity in de- 
tection and parameter estimation. As for detection, the 
SNRs presented here may well be smaller than SNRs for 
PN-corrected eccentric binary inspirals, since the addi- 
tion of harmonics to the amplitudes generally increases 
the power in the signal. As for parameter estimation, 
the eccentricity corrections to the phase of the Fourier 
transform have a different frequency dependence relative 
to PN corrections to the phase in the circular case. This 
suggests that the initial eccentricity might be weakly cor- 
related to other intrinsic parameters. A more detailed 
study is necessary to verify this conjecture, and it will 
be a topic for future work. In Section IVII we outline a 
possible extension of the formalism to higher PN orders. 
When this extension is achieved, ready-to-usc Fourier do- 
main gravitational waveforms could be employed in GW 
searches and parameter estimation. 

The remainder of this paper deals with the details of 
the calculations and results presented above. It is or- 
ganized as follows. Section [IT] presents the basics of the 
Kepler problem and establishes the notation used in this 
paper. Section Hill discusses how to model GWs from ec- 
centric binary inspirals in the time domain, while Sec. II VI 
describes its frequency-domain representation. SectionlVl 
presents the SNR calculation, while Section IVll discusses 



PN corrections. Section IVlII concludes and points to fu- 
ture research. 

Technical details are discussed in the Appendices. Ap- 
pendix 1^ reviews astrophysical scenarios that could pro- 
duce eccentric binaries in the LISA band. Appendices [B] 
and O list some lengthy coefficients appearing in our ana- 
lytic calculations. Appendix ID] discusses the effect of pos- 
sible eccentricity-induced modifications to the innermost- 
stable circular orbit (ISCO), concluding that they are 
negligible in our context. Finally, Appendix [El shows how 
to compute the orbital frequency at some given time be- 
fore merger for eccentric binaries. 

In this work we follow the conventions of Misner, 
Thornc and Wheeler [i^: the metric has signature 
(— , -|-, -I-); spacetime indices are labeled with Greek 
letters, while spatial indices are labeled with Latin let- 
ters; unless otherwise specified, we use geometrical units, 
where G = c = 1, G stands for Newton's gravitational 
constant and c for the speed of light. 



II. THE BASICS OF THE KEPLER PROBLEM 

In this section, we review some of the basic concepts 
related to the Kepler problem in Newtonian mechanics, 
as they are relevant to this paper. We present here only a 
minimal description of this problem and refer the reader 
to for a more detailed account. We also establish the 
notation we shall employ in the remainder of this paper. 

Consider a system of two point particles in an eccentric 
orbit. In the Newtonian Keplerian representation, the 
Newtonian orbital trajectories are given by 
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a (1 — e cos u) , 
I = u — e sin u, 



2 arctan 
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(2.1) 
(2.2) 

(2.3) 



where the notation, following [5l|, is as follows: (p is the 
orbital phase; r is the relative separation vector between 
the compact objects, namely r = r(cos , sin , 0); a is 
the semi-major axis of the ellipse; e is the eccentricity 
parameter; u is the eccentric anomaly; / is the mean 
anomaly; v is the true anomaly, and N is the mean mo- 
tion. The quantities io and 0o are some initial time and 
initial orbital phase that arise as constants of integration. 
Since the energy and angular momentum fluxes depend 
on a and e, and together cause the latter to vary with 
time, it can be shown that a and e co-evolve according 
to @,[| 



- ( J^\'^ - 



a(e) 



cocr(e) , 



(2.4) 



where M = mi -I- m2 is the total mass. The quantity 
F is the Keplerian mean orbital frequency, which can 
be associated with an instantaneous mean orbital fre- 
quency whose evolution is discuss further in Sec. IIV 21 
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The quantity cq is a constant defined by -F(eo) 
the function cr(e) is given by 
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The quantity eg is hencefortli always defined to be the 
eccentricity when the GW signal enters the detector sen- 
sitivity band. For AdvLIGO, ET and LISA, this corre- 
sponds to the initial eccentricity at 20 Hz, 1 Hz or 10 Hz, 
and 10~^ Hz respectively. 

Gravitational waveforms for binary inspirals depend 
on trigonometric functions of the orbital phase, but for 
eccentric inspirals this phase is a complicated function of 
the orbital frequency. In the circular orbit limit (e ^ 0), 
the orbital phase satisfies = 2tt F{t){t — to)j where F 
is the Keplerian orbital frequency (one half the domi- 
nant, quadrupole GW frequency). For eccentric inspi- 
rals, however, the phase is related to the arctangent of 
the eccentric anomaly, which is then related in a tran- 
scendental way to the mean motion, and thus, to the fre- 
quency N — 2'kF . We must then find a way to express 
the orbital phase as a function of the mean anomaly I. 

Let us re-express the cosine and sine of the orbital 
phase in terms of the mean anomaly, through the well- 
known Keplerian relations [4l| 
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Equation (|2.ip allows us to rewrite these relations as 
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Moreover, from the Fourier analysis of the Kepler prob- 
lem, one can expand trigonometric functions of the ec- 
centric anomaly as series of Bessel functions of the first 
kind, Jfc. One then finds that 1411 
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Jfc(fce)cosfcZ, (2.11) 
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where primes stand for derivatives with respect to the 
argument, and 



i-iy 
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with r the Gamma function (see e.g. [53|). 

With these relations, we can now express the cosine 
and sine of the orbital phase as a function of the mean 



— Exact Solution 
■ ■ ■ N=3 
-■ N=7 
■- N=10 
■- N=15 




FIG. 3: Plot of the sine of the phase calculated numerically 
(solid black) and expanded in Bessel functions for a system 
with eo = 0.99, following Eq. (|2.14|l . and keeping 3 (dotted 
red), 7 (dashed green), 10 (dot dashed blue) and 15 (dot-dot- 
dashed orange) terms. 



anomaly. Inserting Eqs. (|2.10p and (|2.1ip into Eqs. 
and (12.91) we find that 



coscj) = -e-f - (1 - e^) ^ Jfc(fce)cosfc/, (2.13) 
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sin0 = (l - e^) [Jfc-i(fce) - Jfc+i(fce)] sinfcL 

fe=i 

(2.14) 



Equations (|2.13p and (|2.14p agree with the correspond- 
ing expressions in the Appendix of Ref. [s^ . These re- 
lations allow us to express the gravitational waveforms 
for eccentric inspirals as explicit functions of the orbital 
frequency. 

The number of terms we should keep in the Bessel func- 
tion expansion depends on the accuracy desired relative 
to the exact solution, as well as on the magnitude of 
the eccentricity. Figure [3] plots the numerical solution of 
Eq. (|2.8p for the sine of the orbital phase with an eccen- 
tricity of 0.99, together with the Bessel expansion of the 
solution given in Eq. (|2.14p . where we keep 3 (dotted), 7 
(dashed), 10 (dot-dashed) and 15 (dot-dot-dashed) terms 
in the sum. Observe that even for such large eccentrici- 
ties we only need fewer than 10 terms to reproduce the 
exact solution quite accurately. 

The number of terms needed to reproduce the exact 
solution does depend on the eccentricity one is trying to 
model. As an example, consider solving Eq. (j2.2p for 
the sine of the eccentric anomaly both numerically and 
with Bessel functions. The latter solution is simply given 
by 111 
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domain. We shall employ the quadrupole formalism, sim- 
ilar to that presented by Moreno-Garrido, Buitrago and 
Mediavilla [3^, \^ , but improved through the techniques 
introduced by Pierro and Pinto [4l|, \^ . 

The starting point is the expression for plus- and cross- 
polarized gravitational waveforms, and hx- We place 
the GW detector at luminosity distance Dl from the 
source, in a direction characterized by the polar angles l 
and (3, defined as those subtended by the local Cartesian 
reference frame of the source and the line of sight vec- 
tor [131 ■ Following Martel and Poisson [s^l, we rewrite 
the expressions of Wahlquist 0| as follows 



FIG. 4; Plot of the absolute value of the fractional relative 
difference between the sine of the eccentric anomaly, calcu- 
lated numerically and expanded in Bessel functions with 9 
terms, for different eccentricities (e = 0.4 solid black curve, 
e = 0.2 blue dashed curve, e = 0.1 red dotted curve). The 
results are here normalized by the exact numerical solution. 



In Fig. m we plot the absolute value of the fractional rel- 
ative difference between the numerical solution and the 
Bessel-expanded solution, keeping 9 terms in the sum. 
We plot the absolute value of the difference, normal- 
ized by the numerical solution for different eccentricities: 
e = 0.1 in dotted red, e = 0.2 in dashed blue, and e = 0.4 
in solid black. Observe that the error due to neglecting 
terms beyond the 9th in the sum amount to less than 
0T% in the worst case (corresponding to the highest ec- 
centricity e = 0.4). For cases with smaller eccentricity, 
such as e = 0.2 and e = 0.1, the relative fractional error 
is always much smaller (10~^ and 10~* for the examples 
above, respectively). 

The small-eccentricity assumption could be removed 
by working directly with the full series in Eqs. (|2.13p 
and (|2T4p . or by resumming it. In fact, Pierro and 
Pinto [4]| have shown how to sum the infinite Bessel 
series to measure the so-called "total harmonic distor- 
tion " which is loosely related to Apostolatos' fitting fac- 
tor [53|. We shall not, however, work with such rcsum- 
mations here, since we shall be interested in binaries with 
eccentricities e < 0.4. In view of this, we shall truncate 
all expressions at order ten in the Bessel expansions. 



III. POST-CIRCULAR EXPANSION OF 
TIME-DOMAIN ECCENTRIC INSPIRAL 
WAVEFORMS 

In this section we describe how to model gravitational 
radiation from eccentric inspiraling binaries in the time 
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- cos(30 - 2/3) + cos(2/?) (1 4- cos t 



/ix = - 



[e cos (j) + e^] sin^ t} , 

^ [4sin(2(/)-2/3)-f-5esin((/)-2/?) 



pDl 

+ e sin(30 - 2P) - 2e'^ sin(2/?)] 



(3.1) 



where = 11111112/ M is the reduced mass, is the lumi- 
nosity distance, (f> is the same orbital phase presented in 
Eq. (|2.3p and p is the semi-latus rectum, which is related 
to the orbital frequency via Kepler's second law 



- = 2^M-i/2 
F 



3/2 



(3.2) 



GWs are clearly dominated, at least for small eccen- 
tricities, by components oscillating at once, twice and 
three times the orbital frequency. Equation (|3.ip is 
valid in the quadrupole approximation, which means 
that higher multipoles (the octupole, hexadecapole and 
higher) have been neglected. These multipoles are pro- 
portional to terms of Oirjc) and higher, which implies 
that Eq. (|3.ip is a good approximation for slow velocities 
and weak gravity. 

The harmonic structure discussed above can be seen 
more clearly if we re-express Eq. p.ip , by using trigono- 
metric identities, as: 
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A 



5e 



C2/3 (1 + C,2) 



5e 



(1 + C?) 



i20 [2c2/3 (1 + c^)] 



sin2(^ [2s2/3 (1 + c- )] + cos 30 



:C2/3 (1 



sin 30 



2 2 I 2/1 

e s, + e 1 



1 -e2 



{cos(/) [— 5es2/3Ci] + sin0 [5ec2f3Ci\ + cos 20 [— 4s2/3Ci] + sin 20 [4c2^Ci 



+ cos 30 [-es2/3Ci] + sin 30 [ec2/3Cj] - 2e2s2;3Cj; } , 



(3.3) 



(3.4) 



r 



where we defined Ci = cos l, Si = sin i, C2^ = cos 2/3 and 
S2/3 = sin 2/3, and we introduced the amphtude 



A 



(3.5) 



with the chirp mass given by = [T'I^M'^I'^ . In the 
hmit e ^ 1 the dominant term is the second harmonic, 
foUowed by the first and third harmonics, while the con- 
stant term contributes only to higher order. 

Explicit expressions for the waveforms as functions of 
time are needed to construct their Fourier transform. 
We shall thus substitute the expansions of the sines 
and cosines of the phase in terms of Bessel functions 
[Eqs. ([2+3]) and ([214]) ] into Eqs. Q and jS^l). The 
Bessel functions, however, are themselves polynomials 
in the eccentricity, which we are also expanding about. 
Three different expansions are thus taking place: 

• PN and Multipole Expansion: Weak-field ex- 
pansion of the metric in terms of mass and cur- 
rent multipole moments of the source distribution, 
which are then expanded in small- velocities. 

• Bessel Expansion: Expansion of the orbital 
phase in Bessel coefficients. 

• Eccentricity Expansion: expansion of Bessel co- 
efficients in small eccentricities. 

The multipolar expansion is a weak-field expansion to 
solutions to the Einstein equations, where we shall here 
keep only the mass quadrupole. This implies that our 
waveforms are accurate only to Newtonian order, where 
we neglect terms of relative order 0{f/c). Such terms 
can be accounted for through a PN analysis, as we shall 
discuss in Sec. IVll 

The Bessel and eccentricity expansions are related, and 
one must ensure that they are performed to a consis- 
tent order. Bessel functions of the first kind behave as 
Jfc(fce) ^ asymptotically for e ^ 1. The phases in 
Eqs. (|2.13p and (|2.14p . however, scale as sin0 ^ cos0 ^ 
Jk{ke)/e ~ e''~^. Thus, an expansion of the waveforms 
to O(e^) requires the phases in Eqs. ((2?T3l) and ((2?T4)) to 
be summed up to fc,„ax = + 1. We shall here work to 
N = 8, which means that the Bessel sums in Eqs. (j2.13p 
and (|2.14p must be performed up to fcniax — 9. 



With these expansions, the waveforms can be written 
as a sum over harmonics. Using that the mean anomaly 
I = 271 F{t — to), the waveforms become 



h+.x = 



where the 



A [c|^.^^ cos (£1) + 5'^^.^^ sin (a) 



e=i 



1 coefficients are 



1 



1 5 1 

e 

192 9216 



(3.6) 



(3.7) 



S_ 



(1) 



(1 + Cj) 



11 

7680 ' 



3 23 o 

r+24^ 



c 



(1) 



S2i3Ci ( 3e - -e 



37 
384* 



19 

256^ 
11 



371 7 
-e 



5120 



3840 



C2/3C 



23 3 



19 , 371 . 

e H e' 

128 2560 



(3.8) 
(3.9) 

(3.10) 



and higher-order coefficients are listed in Appendix |B] 
This expansion resembles that of [40| . but it differs in 
that we are here allowing for arbitrary binary inclinations 
via the angles (i, /3) and we are consistently expanding 
to the same order in eccentricity, without keeping pref- 
actors of (1 — e)^^. We have checked that our results in 
Eqs. (3.7)-(3.10) and those in Appendix[B] are consistent 
with Eqs. (16)-(18) of Ref. [40], after identifying their 
variable O with our l and re-expanding their expressions 
in powers of e. 

The maximum k that one employs in the sums of the 
Bessel expansion, /cmax, is not generically equal to the 
maximum i one uses in the sums of the harmonic decom- 
position of the waveform, £niax- This is because higher 
harmonics in the waveform, such as ft.+,x ^ cos 30, can 
be re-expanded as powers of cos and sin with stan- 
dard trigonometric identities, which then become pow- 
ers of Bessel series via Eqs. (|2.13p and (|2.14p . Cross- 
terms in the product of Bessel series combine to produce 
harmonics of higher order. For example, if fcmax = 4, 
the harmonic decomposition contains terms of the form 
/i+.x ^ /ige^coslOZ. Increasing fcmax leads to terms 
that modify /iq, but only up to fcmax = 9, beyond 
which /iQ is not modified. Generically, this means that 
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1 = TV + 2. so with our choices {N = 8, 
^max = 9), we have Cax = 10. 

The expressions presented above depend on the eccen- 
tricity, which itself is a function of time. Equation (|2.4p 
can be solved for the orbital frequency as a function of 
the eccentricity to obtain F/Fq = [(7{eo)/a{e)]^^^ , where 
_Fo is defined such that F{eo) = Fq. This equation is not 
invertiblc for large eccentricities, but in the limit e ^ 1 
it yields: 



eoX 



-19/18 



1 



3323 
1824^ 



19/9 



+ 



66253974 
" 15994231^ 
105734339801 g 
36410425344 ^° 
2505196889835 
105734339801 ^ 



-19/9 



i-x 

50259743 



15994231 4 
6653952 ^° 

-38/9 



15994231' 
1138825333323 



-X 



-19/9 



-38/9 



105734339801 

1472105896313 
~ 105734339801 ^ 



-19/3 



(3.11) 



where we have defined x = F/Fq. Notice that Eq. (|3ll|) 
is a series in odd powers of eo, and as such, it possesses 
uncontrolled remainders of ©(eg). 



The inversion of the eccentricity as a function of fre- 
quency is not valid for all frequencies and all initial ec- 
centricities. Figure [5] plots Eq. (|3.1ip as a function of 
frequencies for different initial eccentricities, and initial 
frequency Fq = 20 Hz. Observe that for initial eccentric- 
ities eo < 0.6, e{F) decays monotonically as a function 
of frequency, as expected, but for cq > 0.6 it ceases to 
be monotonic, displaying two peaks. This unphysical be- 
havior is a signal that the expansion in Eq. (|3.11[) breaks 
down. This occurs when the first correction in the eg 
expansion ceases to be much less than unity. This re- 
quirement translates roughly to F 3> 1.2 c^^^^Fq, or 
simply eo ^ 0.8. We see then that truncating all expres- 
sions at e < 0.4 is consistent with this requirement. More 
importantly. Fig. [5] shows that Eq. (|3.1ip is well-behaved 
for this range of eccentricities. 

Monotonicity, however, is not sufficient to guarantee 
that Eq. (|3.1ip is valid up to e = 0.4. The ultimate 
test is to compare this power series solution to the ex- 
act numerical inversion of Eq. (|2.4[) . We find that the 
exact numerical solution can be fitted by the following 
phenomenological fraction to better than 1% accuracy: 




200 300 
F[Hz] 



FIG. 5: Eccentricity as a function of frequency for differ- 
ent value of the initial eccentricity eo evaluated at Fq — 
20 Hz. Solid lines correspond to the eccentricity as given 
by Eq. pTTj) for eo = {0.1,0.3,0.5,0.6,0.7,0.8,0.9}, in as- 
cending order. Dotted lines correspond to the approximant 
of Eq. (|3.12p for the same initial eccentricities. 



16.83-3.814/3 



'0.3858 



16.04 + 8.1 /3 



1.637 



(3.12) 



where we have defined P = y^l'^ ja^. One can check that 
in the limit % — > 1 . Eq. (|3.12p equals eo with an accuracy 
of roughly 1%. Equation (|3.12p is plotted in Fig. [5] with 
dotted lines. Observe that the dotted lines are always 
close to the solid lines for eo < 0.7. The use of Eq. (|3.12p . 
however, would go against the philosophy of this paper: 
to power-series expand all quantities in the limit eo <C 1 . 
We leave exploration of this type of resummation and 
other for future work. 



The waveform coefficients can now be written entirely 
as a generalized power-series expansion in the frequency. 
Inserting Eq. (|3.1ip into Eq. (j3.6p and re-expanding in 
the limit e <C 1, wc find, for example 
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(2) 



234273299 



2-5 



2 



19 



2079360 ^ 
1559384621213 

2275651584 
8690279 
415872 ^ ' y 



' 16615 
^ 912"^ 
■ 41434504475 
568912896 ^ 



19 
9 



76 



I eo' 



9 19 

eo X ^ 



19123 , 

469672525907 
^ 758550528 

19 1209 

9 



8448925 



_ 19 
9 



X ^ 



/3323 
1^ 912 ^ 



207936 

778490172577 
632125440 
4 , / 1689785 



63545729 
415872 



19 

X ^ 



(3.13) 



)9 _48\ 4 / 1689785 _i9 1339169 _38 
.A^ J^" ^ V 207936 ^ ' 46208 ^ 
8286900895 _i9 9897925427 _38 28877797117 _u, 1428551432057 _76\ g' 
V 568912896 ^ ' ~ 84283392 ^ ' ^ 126425088 ^ ' " 11378257920 / ^° . 



Notice that in the Umit eo — ^ 
the appropriate circular hmit: 



0, Eq. (|XTI)) reduces to 



C4 



2(l + <)c2/3. As 



we can see, the modified coeflBcients are comphcated and 
uniUuminating, which is why we do not present the re- 
maining ones here. Nonetheless, it is straightforward to 
insert Eq. (|3.1ip into the waveforms of Eq. p.6p to ob- 
tain amplitude corrections as a function of the orbital 
frequency. 



as y ^ -1-00, where '^'■^^(a) stands for the p-th deriva- 
tive of ip evaluated at t = a. In Eq. (|4.2p the fac- 
tor of e^*'^/'-^^-' has a positive sign if ip'^P^a) > 0, and 
a negative sign if il)'^P\a) < 0. Here we have chosen 
the stationary point to be located at i = a, such that 
TA'^^a) = ••• ^^^P~^\a) =0. 

Let us then define the Fourier transform of some time- 
series B{t) as 



IV. FOURIER TRANSFORM OF THE 
WAVEFORM IN THE SPA 

In this section we calculate the Fourier transform of 
the waveform computed in the previous section. To 
do so, we shall employ the SPA (see [1^, and 
for a discussion in the context of GW data analysis), 
which is an expansion in the ratio of the radiation- 
reaction time scale to the orbital period. In this asymp- 
totic expansion we need only keep the controlling factor, 
since subdominant terms can in general be neglected for 
matched-filtering purposes Recently, it has been 

proposed that amplitude corrections in the waveforms 
might play a critical role in the data analysis prob- 
lem Hi, M, m, m, \M, M, M, m, m, [n, iii, but we 

defer a discussion of those corrections to future work. 
We shall here primarily follow the prescription of [38|. 

Let us begin by reviewing the SPA, following Ref. [Sa ]. 
Consider the generalized Fourier integral 



Hy) = / 5(0e^^'^(*)dt, 



(4.1) 



where g{t), '0(i)) b and y are all real. In order to find 
the asymptotic behavior of such an integral as y — > -l-oo, 
one searches for stationary points, namely those where 
-0 = 0. This is because in the neighborhood of station- 
ary points the integrand oscillates less rapidly, and there 
is less cancellation between adjacent subintervals [ssj . 
Thus, the asymptotic behavior of Eq. (|4.ip is given by 



I{y) - gia) f^v^^'^)^^^ li-^v) 



2/|V'(f)(a)| 



i/p 



r(i/p) 



(4.2) 



B{t)e^'"f*dt, 



(4.3) 



and let us write the time-domain waveform as the prod- 
uct of a slowly-varying amplitude A{t) and a rapidly- 
varying cosine with phase i4>{t) and £ > 0. Then, the 
Fourier transform of the cosine (denoted by a subscript 
C) becomes 



Bcif) 



At) (e 



,2Trift+U<p{t) I 2TTift-U4>(t) 



dt. 

. . ^^-^^ 

The first term in Eq. (|4.4p does not contain any sta- 
tionary points and, thus, it vanishes via the Ricmann- 
Lebesgue lemma j55| . The second term, however, 
does have a stationary point at the value to where 
£(t>{to) = 27r/, which defines the stationary phase con- 
dition: F{to) = f/£, where we have defined F{t) = 
(t>{t) / {2tt) . Thus, the asymptotic behavior of the Fourier 
transform of a cosine timc-serics is 



]3c{f) = ■^(^0) p-»(^+7r/4) 



2^j£F{ta) 

where we have defined the phase 



-27r/<o + ^0(io)- 



(4.5) 



(4.6) 



The Fourier transform of a sine times-series is then sim- 
ply Bs{f) = iBcif)- Equation (|4.5p is identical to 
Eq. (I12D with p = 2, g{t) = A{t)/2, y / and 
ip{t) = 2nt — i(j>{t)/f. In obtaining this solution we 
have implicitly assumed that d{lnA)/dt ^ d4>/dt and 
d^(f)/dt^ ^ {d(f)/dt) , which mathematically enforces the 
physical condition that the amplitude varies much more 
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slowly than the phase. In Eq. (|4.5p there is an extra fac- 
tor of two relative to Eq. (|4.2p . because the phase 5" is 
not monotonic in its rangc^. 

In order to find the full solution, we must solve for the 
phase ^I*. Defining the quantity r = F/ F, we can rewrite 
<j){F) and t{F) as 

0(F) = 2^ t' dF\ t{F) = ^dF', (4.7) 



which then leads to 



rF{to) 
nF{to)] = 27r / 



F' 



- dF'. (4.8) 



Of course, these expressions must be evaluated at the 
stationary point, given above by F(io) = f /(■ 



1. Circular Case 

The above formalism can be understood better by 
studying the well-known circular case. Let us then solve 
explicitly for the Fourier transform of the response func- 
tion h{t) in the SPA for a binary in circular orbit (e = 0). 
The response function is defined via the linear combina- 
tion 

h{t) = F+{0s , 05 , i^s)h+ + Fx (Os , <f>s , V'5)/ix , (4.9) 

where F+ and Fx are the so-called beam-pattern func- 
tions that characterize the response of the detector to an 
impinging GW and are also slowly- varying (see e.g. [66l|). 
The orbital frequency evolution is given by 



dF 
'dt 



48 



5ttM'' 



{2ttMF) 



11/3 



(4.10) 



(see e.g. the leading-order contribution to Eq. (A. 2) in 
(4^). We can rewrite the response function of Eq. (|4.9p 
as h{t) = hc{t) + hs{t), where 



The Fourier transform can then be computed in the 
SPA via Eq. gSl) with i = 2. We thus obtain 



1/2 



Dr 



QcicP) [2F{to)] 



-7/6 



X e 



-i(*-|-7r/4) 



(4.14) 



We can solve for the time and phase functions to obtain 

rF 1 



0(F) = (j)c + 2TTj T'dF' 



32 



(27r7V/(F) 



-5/3 



/ t' ^ A// 

;p7'^^' = ic-^(27rXF) 



-8/3 



(4.15) 



where 0^ a-nd tc a-re the orbital phase and time of co- 
alescence. Substituting the stationary phase condition 
F(to) = //2 into these expressions, the phase ^ becomes 



-2^/<, 



128a; 



(4.16) 



where we have defined x = (nAif)^/^, and 0c is the GW 
phase at coalescence. The argument of the exponential 
of the Fourier transform is then 



i In 



'^(circ) 
l^(circ)l 



27: ft, - 0, - J + -L (^^Mf) 



and the full Fourier transform becomes 



-5/3 
(4.17) 



''(circ) 



384 



1/2 



-2/3 



X5/6 

Dl 



Q{i,P)r"' (4.18) 



X exp 



I 2TTft, 



128 



{TTMf) 



-5/3 



where Q — Qc + iQs- Equation (|4.17p is in agreement 
with well-known results in the literature [56| . when we 
keep in mind that the GW frequency / = Few = 2F is 
usually adopted to write down all results in calculations 
involving circular binaries. 



hcit) = AQcii.P) cos 20, (4.11) 
hs{t) = AQs{i,P) sin 20, (4.12) 



2. Eccentric Case 



and where the amplitude ^ is a function of frequency, 
defined in Eq. (|3.5p . We have here introduced the follow- 
ing functions of the polarization and inclination angles: 



QcihP) = 2(l + c2)c2/3F+-4QS2/3Fx,(4.13a) 
Qs{l,P) = 2(l + c2)s2/3F++4c,C2,3Fx.(4.13b) 



Let us now focus on eccentric inspirals. Once more, 
we must consider the response function of the detector, 
defined by Eq. (|4.9p . except that now /i+.x correspond 
to the eccentric waveforms discussed in Sec. IIIIl For ec- 
centric waveforms, the response function becomes 

10 

h{t) ^ A^[Ti cos {£ I) + sin {£ I)], (4.19) 



^ The stationary phase integral can be broken down into two parts, 
inside each of which 5/1 is monotonic. Each of these integrals has 
a stationary phase contribution that leads to the factor of two in 
Eq. 114.21 1. See [S^l for more details. 



where ^ is a function of frequency defined in Eq. p.Sp . 
and where we have defined 



Pf = F+Cf +FxC'l'\ 



E, = F+5f +Fx4'^. 



(4.20) 
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These coefficients are slowly-varying functions of time, 
which can be written as functions of the orbital frequency 
for some given initial eccentricity eq . By using the trivial 
trigonometric identity cos(£ I + (p) = cos(£ I) cos(0) — 
sin(£ I) sin(0) we can combine terms into a single sum of 
the form: 
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Ui cos {i I + (/)£), 



where 



ai ^ sign(ri.)yi^ 



tan M - — 



(4.21) 

(4.22a) 
(4.22b) 



We shall not present the coefficients and 0£ here, but 
they can be straightforwardly calculated using results 
from the previous section. 

The Fourier transform in the SPA then becomes 



384/ Dl 



■2F{to)] 



-7/6 



(4.23) 
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where we have used the fact that for eccentric orbits the 
orbital phase evolution is given by Q 



dF 
'dt 



dadF 
dt da 



{2ttMF) 



11/3 



1 



Z3p2 , 37 4\ 

24 ' 96 / 



o2^7/2 



SttA^^ ' ' (1 -e-, 

(4.24) 

The frequency that appears in Eq. (|4.23p can be thought 
of as an instantaneous mean orbital frequency. Such a 
quantity is "instantaneous" in the sense that it evolves on 
a radiation-reaction timescale. The factor of comes 
about due to the factor of {d^(j)/dt'^)-^/^ in Eq. 
Also note that now ^(to) = ^^(^o) and the factor of 4>e 
cannot be pulled out of the sum because it depends on 
£. One can check that in the limit cq — > 0, a2e~"^^ — > Q 
and we recover the circular limit. 

The phase ^/ must be evaluated at the stationary point 
to, which is here defined implicitly via £l{to) = 27: f or 
simply F(io) = //^, as already discussed. The phase 
is then essentially Eq. (|4.8p , which requires knowledge of 
the characteristic time scale t. Unlike the circular case, 
for eccentric inspirals the integral over r can only be done 
approximately, since 



1 



,7/2 



1 + ie^ 



3Zp4- 

96 



(4.25) 



and e is a slowly varying function of F that cannot be 
inverted in closed form. We can achieve this inversion 



asymptotically for small eccentricities, by first expanding 
Eq. (ji:^ in e < 1 
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(4.26) 



Since the eccentricity as a function of frequency is given 
by Eq. p. lip , the characteristic time becomes 



1 



157 



-X 



1044553 

21888 
3471049619^ 

9980928 ' 



-38/9 



24 
521711 



-X 



-19/9 



-38/9 



21888 

265296245 



-X 



-19/9 



-X 



-19/3 



135641025 

369664 
25654857812777 



4990464 
450735126075 
112377856 



-X 



-19/3 



-X 



18205212672 
1301043440515 .^g/g 
13653909504 ^ 
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54615638016 ^ 
0(ej°)l, (4.27) 



where as usual y_= F/Fq. This is a generalization of 
Eq. (A8) of Ref. [3^ to higher powers of eccentricity. 

We can now compute the new phase [Eq. (|4.8p ] by in- 
tegrating the characteristic time. Using the stationary 
phase condition, F(io) = f/i, we obtain 
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(4.28) 



where we recall that x = (ttA^/)^/"^, and 4* has now be- 
come a function of We have checked that the first few 
terms in the phase of Eg. (I4.28P agree with the phase 
computed in Eq. (AlO) of [33|. Notice that when we ap- 
ply the stationary phase condition to e{F) we must also 
rescale Fg — > /o/^, so that e(/o) = eo. Otherwise, the 
eccentricity function would not be properly normalized. 

Combining all pieces together we obtain the Fourier 
transform in the SPA, namely 
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(4.29) 



12 



where we have defined 
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(4.30) 

-i<l>e(f/k)^ (4 3]^-) 



This is the Fourier transform of the waveform for eccen- 
tric inspirals in the SPA. Note that we have kept up to 
ten harmonics, which corresponds to a consistent expan- 
sion in the eccentricity to O(e^) both in the amplitude 
and in the phase. We aheady saw in Sec. [TTl that this 
is enough to model the Bessel function to high accuracy 
even for relatively high eccentricities. 

The Fourier transform presented here depends on the 
coefficients ^£ that need to be re-expanded in the limit 
Co <C 1. These coefficients can be obtained from 
Eq. (|4.3ip . using the definition of e(F) in Eq. (|3.1ip . ai in 
Eq. (|4.22ap . in Eq. (|4.22bp and Te and S^. in Eq. 
where C+^x and S'+,x a-i'c given in Appendix [Bl The re- 
sulting expression must then be re-expanded in the limit 
eo ^ 1 to ©(cq). We shall not present these expres- 
sions here in full generality, since they are lengthy and 
complicated. Instead we present partial results for as 
a function of e{F) in Appendix [Cl for an optimally ori- 
ented binary (/, = /3 = 0). In the next section, we shall 
employ these expressions in combination with Eq. (j3.1ip . 
and re-expand them in cq «C 1 to eighth order to compute 
the SNR. 



V. SNR CALCULATION 

In this section we compute the SNR using the Fourier 
transform of the waveform in the SPA [Eq. (|4.29p ]. The 
SNR is defined via 
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'-ml „df, 



/ic 
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(5.1) 



where Snif) the one-sided noise power spectral den- 
sity and the star superscript stands for complex con- 
jugation. The noise curves of the AdvLIGO, ET and 
LISA detectors are taken from Refs. [G^I, [111 and (G^I 
respectively; for LISA, in particular, we adopt the sim- 
ple "angle-averaged" model discussed in [66i] . 



A. Limits of Integration 

The upper frequency of integration, /high, is either the 
frequency at which the motion transitions from inspiral 
to plunge or the maximum frequency at which the detec- 
tor noise is under control. Since the noise power spectral 
densities for LIGO and ET increase steeply at high fre- 
quency, we will choose 
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high 
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min [2i^isco, 1 Hz] 
/liigh = 2i^isco- 



(5.2) 
(5.3) 



In the previous equations, as customary in the GW lit- 
erature, we (somewhat arbitrarily) pick the ISCO fre- 
quency to be i^isco = 6~"^/^(27ril/)~^, in analogy with 
the orbital frequency of a test particle at the ISCO of 
the Schwarzschild spacetime. In Appendix ID] we discuss 
possible eccentricity-induced modifications to this con- 
ventional ISCO frequency, concluding that such modi- 
fications should not introduce significant corrections to 
our SNR calculations. 

The lower limit of integration /low is determined by a 
seismic (or acceleration) noise cut-off: 
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(5.4) 



where we shall investigate the ET SNR with /f ^ = 1 Hz 
or ff'^ = 10 Hz. The quantity /acc corresponds to the 
minimum frequency at which acceleration noise is under 
control in LISA. 

Different harmonic components will generically sample 
different frequency ranges if we terminate all integrations 
when the dominant quadrupole GW frequency equals the 
ISCO frequency. To ensure that higher harmonics do 
not exceed the region of validity, following , we shall 
truncate the waveforms with unit step functions Q{x) 
(0(a;) = 1 if a; > and zero otherwise): 
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where we have removed the factor of 7r/4 in the phase, 
since it cancels out in SNR calculations. The step func- 
tion guarantees that higher harmonics are truncated at 
the correct upper frequency cut-off. 

LISA sources can spend several years in the LISA band, 
an issue that must be accounted for, since the detector 
will not take data for more than a few years. Follow- 
ing [3, [6^ , we shall multiply the waveform by an addi- 
tional step function: 



10 



2/3 



-i'i't 



where /yr is the GW frequency of the fundamental har- 
monic at a time T before the system reaches the ISCO 
(see Appendix |E] for a discussion of how to calculate this 
quantity for eccentric inspirals). We shall here choose 
T to be equal to one year (hence assuming, somewhat 
optimistically, that we can observe the whole last year 
of inspiral). This step-function cut-off guarantees that 
all harmonics are integrated for no more than one year, 
which is the higher-harmonic generalization of the crite- 
rion used in [66|. Note also that we have multiplied the 
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LISA waveform amplitude by a geometrical correction 
factor of -\/3/2 (see [HI, [6^ for details). 

With these considerations in mind, the SNR is given 
by 



4SR 



hi 



hA h\ 



df, 



(5.7) 



where A stands for any of LIGO, ET or LISA. Cau- 
tion should be exercised in comparing results between 
different detectors. Even for astrophysical systems with 
the same masses, different detectors have different low- 
frequency cut-offs, and the initial eccentricity cq is de- 
fined as the value of e at that frequency. For example, a 
100 Mq system with eg = 0.3 does not correspond to the 
same astrophysical system when we discuss AdvLIGO, 
whose seismic cut-off is 20 Hz, and when we discuss ET, 
whose seismic cut-off is 10 Hz or 1 Hz. 



B. Results 

Figure [H] plots the SNR for an equal-mass system as a 
function of the total binary mass expressed in solar mass 
units. This SNR is computed at the optimal binary ori- 
entation (l ^ (3 = 9s ^ (f>s = i^s = As discussed 
earlier, for each detector the initial eccentricity is com- 
puted at some (somewhat conventional) lowest cut-off 
frequency, which is different for each detector. We recall 
that for AdvLIGO and ET this lower cut-off corresponds 
to the seismic noise "wall" (20 Hz and either 1 or 10 Hz), 
while for LISA we (conservatively) adopt a lower cut-off 
at 10-4 Hz. 

As a generic trend, the mass reach increases by as much 
as a factor of five for the largest initial eccentricities ex- 
plored here. For systems with cq = 0.4 AdvLIGO could 
observe binaries with total mass up to IO^Mq, ET could 
see systems up to 2 x 10"^ or 2 x 10"* (for a 10 or a 
1 Hz low- frequency cut-off, respectively), and LISA could 
see binaries up to 2 x 10^ Mq. This is to be compared to 
circular inspiral mass reaches of approximately 200Mq 
for AdvLIGO, 4OOAf0 or 4 x IO^Mq (for a 10 or a 1 Hz 
low-frequency cut-off) and 4 x lO^M© for LISA. 

Another feature of this figure is that for low-mass sys- 
tems, the circular SNR curve seems to overlap that of 
eccentric binaries. One can understand this by noting 
that low-mass systems merge in the high-frequency band 
of the detector, since the merger frequency is inversely 
proportional to the total mass. In such cases, the binary 
circularizes before merger. Suppose that the binary's ec- 
centricity reduces to cq < 10~^ at some "circularization 
frequency" Fc- One can then compare the number of cy- 
cles the binary spends in {/iow,^c} relative to the num- 
ber of cycles spent in {Fc, -Fhigh}, to find that the latter is 
overwhelmingly large for low-mass systems. Such a fact 
does not imply that circular waveforms are sufficient for 
detection or parameter estimation to extract signals from 
eccentric inspirals of low-mass. The SNRs shown here are 



"optimal," and thus, a much more careful fitting-factor 
study is necessary to determine whether circular tem- 
plates suffice to extract eccentric binary signals. 

For high masses the SNR presents a somewhat oscil- 
latory behavior. These oscillations seem to scale with 
the eccentricity, becoming worse for systems with cq ^ 
0.4. Oscillations are expected, since different harmonics 
could interfere in the SNR integrand and since the step- 
function truncation of the waveforms will introduce os- 
cillations at overtones of the truncation frequencies. We 
discuss these issues in more detail in the next subsection. 



C. Accuracy of the Approximation 

An important issue concerns the accuracy of our post- 
circular approximation. Our approximation is essentially 
an expansion for eg ^ 1, so it should break down as we 
increase the initial eccentricity. On the other hand, if the 
condition cq 1 is verified, a relatively small number of 
harmonics should model the waveform accurately enough 
that we would not lose much in terms of SNR. 

In order to explore this issue, in Fig. [7] we plot the 
absolute value (5p(^max, 10) = |p(^max) — /o(10)|j where 
p(^max) is the SNR computed by keeping ^max terms in 
the harmonic sum of Eq. (|4.29p . In the top two panels we 
consider systems with moderate initial eccentricity (eg = 
0.01 and bq = 0.1), which are probably most relevant 
for several classes of astrophysical GW sources. When 
eg = 0.01, the deviation in SNR relative to the highest- 
order terms we computed (^max = 10) is at most of 0(1) 
or of 0(10^^) when one uses ^max = 2 and ^max = 3, 
respectively. On the other hand, for the cq = 0.1 case, 
a comparable accuracy in SNR requires ^max > 4 and 
^max > 5, respectively. It should not be surprising that a 
smaller number of harmonics is required for systems with 
low eccentricity. Our analysis suggests that summing up 
to ^max — 4 should bc enough for systems with cq < 0.1, 
while for systems with 0.1 < eo < 0.3 one needs £niax > 8. 

More interesting features emerge for larger values of cq 
(bottom panels in Fig. [7]), which we list below: 

1. Different harmonics play a critical role in the 
SNR at different mass ranges. When eo = 
0.3 the SNR difference peaks at approximately 
(200,400,500,600)Afo for 4iax = (2,3,4,5); for 
lower and larger masses, a smaller number of har- 
monics is necessary. 

2. The number of harmonics necessary to cover the 
entire mass range is a function of the initial eccen- 
tricity cq. When eo = 0.3, for example, harmonics 
with £ > 7 are not needed, because Sp{6, 10) < 7 in 
the entire mass range. On the other hand, for the 
eo = 0.4 case, one really needs at least ^max = 9 
to obtain errors Sp{9, 10) < 10 in the whole mass 
range, while for eg = 0.5 the approximation seems 
to break down, unless more harmonics are included. 



14 




FIG. 6: SNR for an equal-mass binary at optimal orientation as a function of total mass in solar mass units for different initial 
eccentricities. The top figures correspond to the AdvLIGO (left) and ET (right) detectors, while results for the LISA detector 
are shown in the bottom panel. 



3. The oscillations visible in the plots are not neces- 
sarily an artifact of the post-circular approxima- 
tion. Indeed, these oscillations are also present in 
the small eccentricity curves [eg = (0.1,0.2,0.3)] 
of Fig. inland in the top panel [eo = (0.01,0.1)] of 
Fig. [71 If these oscillations were an artifact of the 
post-circular approximation, they would vanish in 
the small eccentricity limit, but instead, although 
they decrease in magnitude, they are still present. 

The fact that different harmonics peak at different 
masses can be understood by observing that two com- 
peting effects control the SNR difference: the eccentric- 
ity decay and the frequency band over which we perform 
the integration. For small masses, one is integrating over 
a larger frequency band, and the binary rapidly circular- 
izes before merging. Less harmonics are needed in the 
limit of very small mass, since the binary is essentially 
circular before reaching the most sensitive region of the 
detector. On the other hand, for really high masses one 
is integrating for short times and essentially capturing 
only the behavior near the ISCO. In such cases, the SNR 
difference is converging to zero, because the SNRs them- 



selves are essentially vanishing (i.e. the range of integra- 
tion asymptotes to zero). 

The oscillatory features can be understood by studying 
the analytic structure of the waveforms used to compute 
the SNRs. For small eccentricities, the oscillations are 
probably due to interference between the different har- 
monics and to the use of step functions to truncate the 
SNR at different harmonics. The waveform contains a 
sum of ten different oscillatory functions, and when this 
sum is multiplied by its complex conjugate, one naturally 
obtains interference of the type exp[i{£ — £')t]. More- 
over, the step function truncation of the SNR also forces 
oscillations at overtones of the ISCO frequency. For ex- 
ample, Figures |S1 and [71 show oscillations at £i^iscoj with 
£ = {1, 10}. This is because the leading harmonic has a 
mass reach corresponding to twice the ISCO frequency, 
whereas the i-th harmonic has a mass reach correspond- 
ing to £ times the ISCO frequency. Since we use a step 
fimction cut-off for every harmonic [see e.g. Eq. (|5.5[) ]. 
the resulting SNR will show the "bumps" seen in Fig. [71 
Of course, although there arc analytic reasons that ex- 
plain the presence of these oscillations, one cannot for- 
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mally exclude the possibility that (for large Cq) some of 
these oscillations are induced by inaccuracies in the post- 
circular approximation. 



D. Angular dependence 

In all SNR plots, we have so far assumed that binaries 
are optimally oriented. This is a very special configura- 
tion, and it is important to investigate the variations in 
SNR for non-optimally oriented binaries. Figure [5] shows 
histograms of the SNR distribution under the hypothe- 
sis of a uniform distribution of the angles involved: the 
angles (t, /3) describing the orientation of the source and 
the three angles {9s,(f>s, ips) appearing in the antenna 
pattern functions [6^1 ■ We binned the data of the 1000 
random realizations in bins of 50. The histograms are 
representive of the number of realizations in each bin. 
Since they are normalized, the numbers on the vertical 
axis do not correspond to the actual number of realiza- 
tions in each bin. The comparison of distributions corre- 
sponding to various eccenricities makes sense only if the 
histograms are normalized. In the two panels of Fig. [8] 
we plot the AdvLIGO SNR distribution for two represen- 
tative binaries with total mass IOGM0 (left panel) and 
300 Mq (right panel). The SNR of the 100 Mg system is 
dominated mostly by the leading harmonic, whereas the 
300 Mq system has significant contribution from higher 
harmonics (see the top panel of Fig. [6|) . In each panel we 
show three histograms, corresponding to initial eccentric- 
ities Co = 0.1 (dotted black), 0.2 (dashed blue) and 0.4 
(solid red). 

The effect of eccentricity is to shift the SNR distribu- 
tion to higher masses. For the IOOMq system, although 
the SNR distribution visibly shifts to the right, the ef- 
fect is rather mild for the initial eccentricities considered 
in the plot. In this case, the distributions look more 
or less similar, except for a longer tail at large SNRs 
for large initial eccentricities. For the 3OOM0 system, 
the shift in the distribution is much more pronounced. 
While the distribution is narrowly peaked at low eccen- 
tricities, it becomes much wider for larger eccentricities. 
This widening of the distribution might be because, for 
large eccentricities, the waveforms are more sensitive to 
terms proportional to q = cos/, (see the expressions of 



C4" and in Appendix [B| . These amplitude cor- 
rections, being proportional to the eccentricity, are sup- 
pressed for small values of cq. 

The SNR histograms show that, although the SNR val- 
ues for any particular observation may deviate consider- 
ably from the optimal values we have quoted, the general 
trends in the SNR induced by eccentricity (mass reach in- 
crease, accessible volume increase) will not be modified. 
That is, regardless of the location of the source in the 
sky, eccentricity has the net effect of increasing the SNR 
or mass reach of any given signal. 



VI. PN CORRECTIONS 

Although this paper has only considered Newtonian 
waveforms to exemplify the post-circular construction of 
the SPA Fourier transform of eccentric waveforms, it 
is instructive to discuss how to extend these results to 
higher PN order. In this section we shall mainly follow 
the conventions of Gopakumar and Iyer 51|, where a 
"2PN term" is one of C'(r/c)'* smaller than the leading- 
order term, i.e. we employ a relative order counting 
scheme. 

When eccentric orbits are studied beyond IPN order, 
one discovers that the Keplerian parameterization must 
be corrected. In particular, Eqs. (|2.1[) and (|2.3p arc en- 
hanced by new Oir jcf' corrections, while Eq. (|2.2p is en- 
hanced by corrections of 0{r/c)'^. Moreover, if we want 
to keep a Keplerian-inspired parameterization we must 
introduce three distinct eccentricity parameters: mea- 
sures radial oscillations, measures azimuthal oscilla- 
tions and ct measures the frequency eccentricity. These 
eccentricity parameters can be written as functions of et 
in a PN expansion, but this choice is somewhat arbitrary 
(see e.g. [701 ^r comparisons of the different definitions 
with numerical relativity simulations of eccentric merg- 
ers) . 

The waveforms then acquire two sets of modifications: 
amplitude and phase corrections. The amplitude correc- 
tions, as expected, take the form 
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where C, = mN is a natural PN expansion parameter 
of 0{r/c)'^, 77 = fim is the symmetric mass ratio and 

(n) 

X are products of functions of the eccentricities and 
harmonic trigonometric functions. The phase corrections 
can be split into 



0-(/.o = A(O + M/(O, 



(6.2) 



where A(Z) is a 27r if -periodic function of the mean 
anomaly, while W{1) is periodic in /, and thus, 27r- 
periodic. The quantity K measures the advance of the 
periastron per orbital revolution [g^l, while the effect of 
W{1) is to modulate the amplitude via nutation. 

The leading order correction to the phasing of GWs is 
due to pcriccnter precession. This effect is embodied in 
the function X{1), which can be written as 



A(0 = Kl = [1 + kpiet)] I, 



(6.3) 



where we have defined K ^ 1 + kp{et). To this order, the 
mean anomaly continues to be given by Eq. (|2.2p with 
e — > et, while the precession correction kp{et) is given by 
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Since A'' appears here as a IPN order correction, we 
can take its Newtonian value in this equation, namely 
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FIG. 7: Absolute value of the difference in SNR computed up to l-max = (2, 3, 4, 5, 6, 7, 8, 9) and up to ^max = 10. We use the 
AdvLIGO noise curve and an initial eccentricity of eo = 0.01, eo = 0.1, eo = 0.3 and eo = 0.4 (top to bottom). The insets 
zoom on the i/-axis in the region between 200 and 800 solar masses. 



N = 27r/, such that Q ~ 2itMF. Thus, the precession 
correction becomes 



, , , 3 (27rAfF)^/^ , ,4 
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(6.5) 



We can expand this function for small eccentricities e* ^ 
1 to find 

A(/) - / [l + 3 (1 + + + + e^)l (6.6) 

The second correction to the phase is given by the nu- 
tation fimction W , defined by 

W = [v - u + et&vau)[l + kp{et)] + 0{f / cf , (6.7) 

where the true anomaly is given by 
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Note that the true anomaly depends on and not et, 
but these quantities are related via 



(6.9) 



In Sec. [IT] we already discussed how to solve for m as a 
function of I in terms of a series of Bessel functions. In 
particular, one can show that 



Z + ^ ( - j Js{se) sin(sZ). 



(6.10) 



Equations (|6.9p and (|6.10p can be substituted in Eq. 
to find V as a function of I. Then Eq. (|6.7[) yields W as 
a function of /. Expanding in 64 ^ 1 one finds 



W{1) ~ e4-(-10 + r?)C^/''^ + 2j sin(0 

-1/2 (4?7- 31)C^/3 + 5/2! cos(/)sin(0 

(-1/6 (-186 + 2777) + 13/3) cos^ [1) 

- 1/6 (12-6?7)C^/3-4/3jsin(0 + ... (6.11) 

This function, however, is part of the phase, so it enters 
the waveform as the argument of trigonometric functions. 
Note that W{1) is linear in et, and thus, when the cos(0) 
or sin((/)) are expanded in ^ 1, W{1) introduces higher 
harmonics into the waveforms. 
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FIG. 8: Histograms of the AdvLIGO SNR (see text) for a set of random observer orientations. The left (right) panel corresponds 
to a M = IOOM0 {M = 3OOM0) system. Systems with initial eccentricity eo = (0.1, 0.2, 0.4) are shown in dotted black, dashed 
blue and solid red, respectively. 



We see then that to IPN order, it suffices to consider 
the pericenter precession correction through A. The cor- 
rections produced by W{1) are automatically accounted 
for in the Bessel expansion. In essence, this is because 
Eqs. (|2.1[) and (|2.2p arc not modified to this order. The 
IPN-corrcctcd waveforms are then (schematically) 



h{t) = ai cos I 
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(6.12) 



where the tkS are constants. If we were to consider 2PN 
corrections to the waveforms, then the formalism out- 
lined here would have to be extended and W{1) would 
contribute by introducing new corrections not accounted 
for in the Bessel expansion. 

The structure of the IPN time-domain waveform in 
Eq. ([O^ is different from that obtained in Refs. [li,!!^], 
in that the above equation does not lead to periastron- 
precession side-bands in the GW spectrum. In Refs. [H, 
l40j such side-bands arise due to the assumption that pe- 
riastron precession leads to a constant 7 oc fcp. This 
assumption breaks down on long time-scales as perias- 
tron precession is not constant, an effect one can justly 
treat as a 1 PN contribution. As a result, one loses the 
artificial side-band structure in the GW spectrum. The 
implications of this effect will be assessed in future work. 

PN corrections modify the Fourier transform of the 
waveform in the SPA. The phase is modified by a 
factor oi {\ + kp), and we now obtain 

^,{F) = kX[t{f/e)]-2nft{f/£), (6.13) 



where 
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The A term contains the (1 -I- kp) dependence that we 
referred to via Eq. (j6.3p . All the machinery developed in 
the previous section then carries through, with the proper 
enhancement of the Newtonian waveform to higher PN 
order. The net effect of higher PN corrections in the 
Fourier transform is to introduce an infinite set of har- 
monics and PN corrections to the lower-order (Newto- 
nian) harmonics considered earlier. 

While considering higher PN order effects, one can 
also work with the PN parameter (Mw)^/^, where w is 
the orbital frequency. As pointed out in Ref. [IH, this 
parametrization helps to more easily recover the circular 
limits of various elliptic-orbit expressions. Furthermore, 
in comparing numerical relativity results to PN expan- 
sions, Ref. [Til found that waveforms parametrized in 
terms of (Muj)'^^^ are in better agreement with numeri- 
cal waveforms. This deserves more careful study in the 
future. 



VII. CONCLUSIONS 

We have proposed a new scheme, the post-circular 
approximation, to construct "ready-to-use" , analytic 
Fourier-domain gravitational waveforms produced by ec- 
centric binary inspirals. The scheme consists of expand- 
ing all quantities in a power series about zero initial ec- 
centricity. We find that the first 10 terms in the Bessel 
solution to the Kepler problem suffice to reproduce the 
eccentricity evolution to better than 0.1% for eccentrici- 
ties e < 0.4. The resulting waveforms are then rewritten 
in terms of ten physical parameters (the reduced mass 
77, the initial eccentricity eo and frequency Fq, the to- 
tal mass M, the luminosity distance £>l, four angles i, 
/9, 9, 4> describing the relative orientation of the source 
and detector, and a polarization angle ip) and the orbital 
frequency, which can be thought of as a function of time. 

This scheme allows us to analytically construct the 
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Fourier transform of the response function through the 
SPA, where one assumes that the radiation-reaction time 
scale is much larger than the orbital time scale. The re- 
sulting Fourier-domain waveforms contain eccentricity- 
induced, higher-harmonic amplitude and phase correc- 
tions. By computing the SNR as a function of total mass 
and eccentricity we find that the amplitude corrections 
increase the mass reach of the detectors by a factor ~ 5 
for moderately eccentric systems, which in turn implies 
that the source volume accessible to the detectors would 
be increased by almost two orders of magnitude. 

The results presented here cannot be used directly in 
realistic data analysis pipelines because PN corrections 
to the amplitude and phase have not been included. In- 
stead, the present paper was concerned with proposing 
a method to construct "ready-to-use" , analytic expres- 
sions for the Fourier transform of the response function, 
which was exemplified through Newtonian-accurate ex- 
pressions. Future research should include such PN cor- 
rections. 

Another interesting research direction is the study of 
the effect of eccentricity in parameter estimation. Ec- 
centricity adds more complexity and information to the 
waveforms that could break parameter degeneracies, thus 
possibly leading to better accuracy in parameter estima- 
tion. On the other hand, the inclusion of eccentricity- 
induced corrections to the GW phase could mimic certain 
high-order PN phase corrections, which would then cre- 
ate new degeneracies. A more detailed parameter estima- 
tion study is needed to assess whether GW measurements 
will benefit or not from the inclusion of eccentricity. 

While finalizing the draft we learned that two different 
groups arc investigating frequency-domain gravitational 
waveforms for eccentric binaries (t^ . pTSj . It would be 
interesting to compare their approach with ours and with 
the time-domain waveforms of 1301. 
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APPENDIX A: LISA ECCENTRIC BINARIES 

In this Appendix we briefly review some literature 
on scenarios leading to non-eccentric binary inspirals in 
the LISA band. We consider in turn stellar mass bina- 
ries, extreme- and intermediate-mass ratio inspirals (EM- 
RIs/IMRIs, respectively), and the coalescence of massive 



BHs. 



1. Stellar mass binaries 

Many stellar mass binaries involving neutron stars are 
expected to be eccentric in the LISA band (see e.g. [ssj). 
It is also well known that LISA should provide a large 
observational sample of interacting white-dwarf binaries, 
whose evolution is driven by radiation reaction, tides and 
mass transfer [t^I . It was recently realized that eccentric 
double white dwarfs formed in globular clusters would be 
detectable by LISA out to the Large Magellanic Cloud 
[t^ . In these binaries, the periastron precession has con- 
tributions due to general relativity, but also to tidal and 
rotational distortions. Tides and stellar rotation should 
dominate at frequencies above a few mHz. The Fisher- 
matrix analysis of [76| pointed out the interesting possi- 
bility to study white dwarf structure with LISA. However 
their analysis neglected the contribution of radiation re- 
action effects, that should be relevant for / > 0.5 mHz. 
We expect our post-circular formalism to be useful in 
this context, since radiation reaction in these eccentric 
binaries should be well modeled by the quadrupole ap- 
proximation. 

2. Extreme and intermediate mass ratio inspirals 

Formation scenarios for EMRI and IMRIs, involving 
a SMBH and either a compact stellar mass object or an 
intermediate-mass BH (IMBH) , are reviewed in Ref. [tJ . 
If an EMRI is formed following a tidal binary separation 
event, the compact star is deposited on an orbit with 
semi-major axis ~ 10^ — 10^ AU and e ~ 0.9 — 0.99, 
and the orbit should circularize by the time it enters the 
LISA band. However, typical EMRIs are expected to 
form by scattering of the compact object into nearly ra- 
dial orbits followed by inspiral due to dissipation, and in 
particular due to GW-emission. Hopman and Alexander 
[Tq showed that the eccentricity distribution of EMRIs 
is skewed to high-e values, with a peak at e ~ 0.7, at an 
orbital period of ~ 10^ s. The dynamical evolution of 
IMBH binaries formed in dense stellar clusters, using a 
combination of A^-body simulations and three-body rela- 
tivistic scattering experiments, shows that the eccentric- 
ity of these systems in the LISA band can be as large 
as ~ 0.2 — 0.3 [t^. The post-circular approximation de- 
veloped here could be applied to these systems once PN 
corrections are taken into account. 



3. Massive black hole coalescence 

The eccentricity of SMBH binaries has been the sub- 
ject of some debate. Gravitational radiation reaction 
alone is not sufRcient to produce mergers between mas- 
sive BHs, which probably require dynamical interactions. 
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Analytic calculations and iV-body simulations show that, 
in purely coUisionless spherical backgrounds, the ex- 
pected equilibrium distribution of eccentricities is skewed 
towards high e ~ 0.6 — 0.7, and that dynamical friction 
does not play a major role in modifying such a distri- 
bution (see Ref. [801, particular Fig. 5). The actual 
eccentricity of a merger event is therefore determined 
by the competition between dynamical wandering and 
GW-induced circularization. Reference [tJ presents ar- 
guments supporting circularization of most binaries by 
the time they enter the LISA band. However, several 
mechanisms producing non-zero eccentricity have been 
proposed in the past (see e.g. Section 2 of Ref. [8l|). 

Recent smoothed-particle hydrodynamics simulations 
follow the dynamics of two BHs orbiting in massive, ro- 
tationally supported circumnuclear discs [s^ . [ssj . The 
rotation of the disc circularizes the orbit if the pair coro- 
tates with the disc. Circularization is efficient until the 
BHs bind in a binary, though in the latest stages of the 
simulations (when the separation is of the order of a few 
parsecs) a residual eccentricity e > 0.1 is still present. 
Circularization possibly reduces the gravitational radia- 
tion merging time scale so much that the binary stalls, 
and no coalescence results. For corotating discs, the nu- 
merical resolution of the simulations is not sufficient to 
compute the residual eccentricity when the BHs are close 
enough that gravitational radiation takes over. More- 
over, if the orbit of the pair is counterrotating the initial 
eccentricity does not decrease, and BHs may enter the 
GW-dominated phase with high eccentricity. 

Collisional processes (such as three-body encounters 
with background stars) may become important at BH 
separations < 6 pc, possibly leading to an increase in 
eccentricity balancing the circularization driven by the 
large-scale action of the gaseous and/or stellar disc. Sev- 
eral investigations show that eccentricity evolution may 
still occur in later stages of the binary's life, because 
of close encounters with single stars [SJ] and/or gas- 
dynamical processes [s^. In particular, the gravitational 
interaction of the binary with a surrounding gas disc is 
likely to excite BH binaries to eccentricities e > 0.1. The 
transition between disc-driven and gravitational wave- 
driven inspiral can occur at small enough radii that a 
small but significant eccentricity survives, with typical 
values e ~ 0.02 (and a lower limit of e ~ 0.01) one 
year prior to merger (cf. Fig. 5 of [11]). If the bi- 
nary has an extreme mass ratio q < 0.02 the residual 
eccentricity can be considerably larger (e > 0.1). Re- 
cent simulations by Cuadra et al. |86[ investigate the 
evolution of the orbital parameters of binaries embedded 
within geometrically thin gas disks. For binary masses 



10^ Mq <M< IO^Mq, they find that orbital decay due 
to gas disks may dominate the binary dynamics for sepa- 
rations below a ~ 10^^ — 0.1 pc, and that in the process 
the eccentricity grows at a rate de/dt ~ 1.5 x lO^^Worb, 
where Worb is the orbital frequency. Saturation of the ec- 
centricity growth is not observed up to values e > 0.35, so 
the binary may have significant eccentricity by the time 
gravitational radiation takes over. 

Stellar dynamical hardening might also leave the bi- 
nary with non-zero eccentricity. Early studies suggested 
that any such eccentricity would be small [m, ^ST,, ,8^ 
(but see and [l^l for examples of eccentricity growth 
in iV-body simulations). More recent TV-body simula- 
tions combined with a Fokker-Planck model [9lj find that 
perturbations of the (initially circular) binary orbit from 
passing stars produce significant eccentricity around or 
even before the time when the binary becomes hard. The 
averaged eccentricity growth is maximum for equal-mass 
binaries with e ~ 0.75 and falls to zero at e = and 
e = 1. It is hard to estimate the final eccentricity, which 
strongly depends on noise-induced changes in e at early 
times, and would presumably be much smaller than the 
simulations suggest in the large- regime of real galaxies. 

Berentzen et al. [48j present simulations following the 
SMBH evolution in rotating galactic nuclei from kpc sep- 
arations down to coalescence, including post-Newtonian 
(PN) corrections to the binary equations of motion. They 
find that the orbital eccentricities remain large (between 
0.4 and 0.99, with typical values around e ~ 0.9) un- 
til shortly before coalescence, and that higher harmon- 
ics of the eccentric signal are detectable by LISA with 
large SNR. Most of these binaries have sizable eccentric- 
ities (up to ~ 0.2) by the time they reach a separation 
~ 10^ Schwarzschild radii, which roughly corresponds to 
a binary of mass M « 2 x 10^ Mq entering the LISA 
band (see e.g . Figure 8 of [11]). The study by Scsana 
et al. [93, 193I [9^ confirms that binaries with mass ratio 
q = M2/M1 < 0.1 and/or eccentricity e > 0.3 can shrink 
to the GW-dominated regime within a Hubble time (see 
in particular Figure 7 of [93]; Section 4.1 and Figure 10 
of [9J]). Last but not least, an interesting scenario pro- 
ducing highly eccentric mergers that could be observed 
by LISA involves close triple SMBH encounters [H, . 



APPENDIX B: HIGHER-ORDER COEFFICIENTS 

In this Appendix we list some of the higher-order co- 
efficients appearing in the expansion (|3.6[) . 
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APPENDIX C: THE (,k COEFFICIENTS 



The coefficients defined in Eq. (|4.3ip are given by the 
foUowing expressions, when one fixes /3 = t = 0: 
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APPENDIX D: ENDING FREQUENCY FOR 
ECCENTRIC BINARIES 

In this Appendix wc discuss possible generalizations of 
the notion of an ISCO to eccentric binaries. The idea 
is that eccentric binaries will transition from inspiral to 
plunge at a frequency slightly different from the circular 
ISCO frequency, and we may worry about the effect of 
this modified ISCO on the upper cut-off frequency used 
in SNR calculations. 

A possible way to modify the ISCO location is to use 
the Newtonian formula in Eq. (|2.4[) with a ~ p/{l — e^) 
and p = 6 + 2e, which corresponds to the value of the 
separatrix between stable and unstable (plunging) orbits. 
In this way we would find that the ISCO frequency is 

1 /l-e2\'/' 

This guess cannot be valid for large eccentricities, when 
the pericenter becomes small, since then the Newtonian 
relations break down. A more accurate approximation 
of the ISCO frequency is to use the pericenter frequency 
57p = M/rp3 at the separatrix pericenter Tp = 6 + 2e, 



leading to [13, [H, ^ 

1 / 1 + e 
^■^<=°==2^(,6T^j ' 

but this result is also not appropriate here because 
this eccentricity corresponds to that associated with 
Schwarzschild geodesies, so it is not equivalent to the 
Newtonian definition of eccentricity we have used (see 
e.g. [lO] for a discussion). 

One expects the residual eccentricity any binary could 
have by the time it enters the strong field to be small. 
The classic work of Peters Q and Peters and Mathews 
suggests that a binary with some moderate initial eccen- 
tricity will rapidly circularize. Since e/en ^ (///o)~^^^^^ 
to leading order [see e.g. Eq. (2.34) in [lOl], an orbit with 
initial eccentricity eg = 0.4 at the beginning of the LIGO 
band will have a final eccentricity of e ~ 0.035 by the time 
it reaches LIGO's highest sensitivity region at 200 Hz. 
These results suggest that the ISCO frequency for eccen- 
tric inspirals will generically be close to the ISCO fre- 
quency for circular inspirals, provided this frequency is 
much larger than the initial frequency associated with the 
initial eccentricity. If the latter is not the case (e.g. if a bi- 
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nary with eo = 0.4 at Fq = 20 Hz merges at 40 Hz), then 
one might have to worry about the precise definition of 
the ISCO, but in such cases the SNR will be dominated 
by the merger waveform and not the inspiral. We are 
thus justified to ignore eccentric corrections to the ISCO 
and employ the usual circular-orbit ISCO expression in 
our SNR calculations. 



APPENDIX E: FREQUENCY AT A GIVEN TIME 
BEFORE MERGER 



LISA sources can easily orbit for more than one year 
in the LISA band. The LISA mission, however, is not ex- 
pected to last for more than a fcw years in orbit. For this 
reason, it is customary to perform LISA SNR and param- 
eter estimation calculations assuming that the source is 
observed over the last year (or few years) of inspiral. 

In the case of circular inspirals, one can compute ex- 
actly (to Newtonian order) the frequency at a given time 
T prior to merger. This is given by Eq. (2.15) in Ref. [6^ : 
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This equation can be obtained by finding T{F) as the in- 
tegral of and then inverting the resulting expression 
to find F{T). 

In the case of eccentric inspirals, an analogous rela- 
tion cannot be obtained analytically. This is because the 
equation for F in Eq. (|4.24|) is a function of the eccen- 
tricity, which itself is a function of the frequency (and 
implicitly time). One could attempt to construct an ap- 
proximation for F{T) by inserting Eq. (|3.1ip for e{F) into 
the expression for F in Eq. (|4.24p to compute T{F), and 
then perturbatively inverting this relation to find F{T). 
The resulting asymptotic scries, however, is poorly con- 
vergent for large masses or large integration times. 

A numerical procedure is thus necessary to find /yr for 
eccentric inspirals. One such scheme is as follows. Given 
some Co and Fq, one can find the corresponding initial 
semi- major axis ag from the first equality in Eq. (|2.4p . 
From this, one can then use this same equation to find 
the corresponding constant Cq. The eccentricityCT a time 
T before merger is then given by Eq. (5.14) in [8|, namely 



T(ao,eo) = --y^ 



(1 
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121 
304* 
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where f3 = mim2M . This is because the eccentric- 
ity at zero orbital separation (roughly corresponding to 
"merger") vanishes in the Newtonian approximation. If 
we set T = 1 yr, then we can solve Eq. (|E2[) for Cyr nu- 
merically using bisection or the secant method in Math- 
ematica. Once the appropriate eyj. is found, one can use 
this to find Cyr via Eq. ()2.4|) . which can then also be used 



to find /, 



for some given eo and A/. 




10 
M[MJ 

FIG. 9: Frequency one year prior to merger as a function 
of total mass for equal mass binaries with different initial 
eccentricity: eo = (crosses), eo = 0.1 (circles), eo = 0.2 
(squares) eo = 0.4 (diamonds). The quantity / — 10~* Hz 
corresponds to the acceleration noise cut-off frequency. 



Figure [9] plots the dominant GW frequency (twice the 
orbital frequency) one year prior to merger found with 
the above algorithm as a function of total mass for equal 
mass binaries with different initial eccentricities. Observe 
that as the eccentricity increases, /yr decreases faster 
with total mass than in the circular case, because ec- 
centricity speeds up the inspiral. Thus, given a fixed 
inspiral time (e.g. one year), the starting frequency must 
be pushed to lower values. Care must be taken, however, 
since /yr appears multiplied by £/2 in the step functions 
used to truncate the waveform. Nonetheless, the above 
figure suggests that for total masses M > lO^M© it does 
not matter whether one uses the circular or eccentric ex- 
pressions for /yr. 
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